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ABSTRACT 

Pulsar motions in our Galaxy from their birth until 2 Gyr are studied statistically 
via Monte-Garlo simulation of 2 x 10^ pulsars with the best available representation 
of the Galactic potential. We find that the distribution of height above the Galactic 
plane for pulsars with characteristic ages less than about 8 Myr could be well fitted 
by a Gaussian function. For older pulsars, an extra exponential function is necessary 
to fit the distribution. The scale-height of the Gaussian component increases linearly 
with time until about 40 Myr. The height distribution becomes stabilized after about 
200 Myr. These results are not sensitive to initial height or radial distributions. Taking 
the relationship between the initial velocity and height distribution, we found from 
the latest pulsar catalog that the height distribution of pulsars younger than 1 Myr 
directly implies the mean initial velocity of 280 ±96 km s~^. Comparison of simulated 
sample of pulsars with the current available millisecond pulsars shows that their ID 
initial velocity dispersion should be most probably 60 ± 10 km s^^. 
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1 INTRODUCTION 



Pulsars are high velocity objects in our Galaxy (Lyne & 
Lorimer 1994). They were born in supernova explosions 
near the Galactic plane where their progenitors live, but 
move away very fast from the plane. The evolution of pulsar 
heights was considered (Bhattacharya et al. 1992; Hartman 
et al. 1997; Mukherjee & Kembhavi 1997) but not explicitly 
described in most pulsar population syntheses, therefore, the 
picture for pulsar height evolution was not very clear. We 
want to demonstrate in this paper by 3D simulations how 
pulsars move in our Galaxy. 

Birth velocities of both normal and millisecond pulsars 
(MSPs) have been investigated by many authors. For normal 
pulsars, Lyne & Lorimer (1994) studied transverse speeds of 
29 pulsars younger than 3 Myr and obtained the mean birth 
velocity Vb =450±90 km s~^. Taking into account selection 
effects and using all available pulsar proper motion data 
in population synthesis, Lorimer, Bailes & Harrison (1997) 
found the mean birth velocity Vb ~500 km s"^ . Hansen & 
Phinney (1997) considered the distribution of proper mo- 
tions, including all available upper limits, and concluded 
that the mean birth velocity Vb should be smaller, only 
250-300 km s"\ Cordes & Chernoff (1998) performed a 
detailed analysis for measured velocities of 49 young pul- 
sars. They favoured a two-component Gaussian model in 
3-dimension with characteristic velocities* of 175 km s^^ 



In the Maxwellian distribution, velocity in each dimension is 



and 700 km s~^. Arzoumanian, Chernoff & Cordes (2002) in- 
fer the velocity distribution of radio pulsars detected in 400- 
MHz surveys, taking into account beaming, selection effects 
and luminosity evolution. They also favour a two-component 
Gaussian model with Gv ~ 90 km s^^ and 500 km s^^. 

For MSPs, Lorimer (1995) accounted for known sur- 
vey selection effects and obtained the local surface density 
of MSPs, birthrate and the lower limit of the mean birth 
velocity that is (t„ > 50km s~^, i.e. Vb > 80km s~^. The 
results were confirmed by Cordes & Chernoff (1997), who ob- 
tained the velocity perpendicular to the Galactic plane, Vz = 
52tjjkm s~^, i.e. a 3-D velocity dispersion (t„ = 84 km s~^ 
using likelihood analysis on previous surveys for MSPs plus 
selection effects. Later, Lyne et al. (1998) found that the 
population syntheses of the MSPs with the Maxwellian ini- 
tial ID velocity dispersion of 80 ±20 km s~^ agrees the best 
with the data of newly discovered MSPs in Parkes southern 
sky survey together with previous MSPs. The 3D mean birth 
velocity is correspondingly 130 ± 30 km s~^. These results 
were later confirmed by Toscano et al. (1999) using more 
MSP proper motions obtained from timing measurements. 

Previously, Narayan & Ostriker (1990) tried to model 
the observed pulsar populations. They deduced an analytic 
formula for the height distribution by assuming a Gaussian 
distribution at all times and solving the energy conserva- 



a Gaussian distribution. The mean velocity in 3-D and velocity 
dispersion in 1-D are connected by Vb. ~ o"t, -v/S/tt ~ 1.6(t„. 
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tion equation in 0-direction (perpendicular to the Galactic 
plane). However, the height distribution changes from time 
to time when pulsars run away from the Galactic plane. The 
integrated ^-distribution of pulsars with different ages may 
not be a Gaussian. We will show by our Monte-Carlo sim- 
ulation that the scale-height evolution is much more com- 
plicated than a Gaussian. Analytic solutions obtained from 
the over-simplified form of a Galactic acceleration model can 
hardly represent the realistic situation, especially over long 
time with non-linear evolution of pulsar heights. By simula- 
tions, we are able to track the height distribution of a pulsar 
sample, for any form of acceleration. 

In our work, the initial velocity dispersion for MSPs is 
obtained by comparison of the height distribution between 
simulated and observed pulsar samples. We also tried to de- 
rive the initial velocity dispersion of normal pulsars by com- 
paring their scale-heights and characteristic ages. In Section 
2 we will describe the simulation procedures together with 
input parameters and governing equations. Several possibil- 
ities of parameter choice are also discussed. In Section 3, we 
show the simulation results, mostly regarding how the pul- 
sar scale-height evolves. We demonstrate in Section 4 two 
applications of our simulations, i.e. to determine the birth 
velocity dispersions of normal pulsars and MSPs. The influ- 
ence of the selection effects is also discussed. The conclusions 
are given in Section 5. 



2 DETAILS OF SIMULATION 

In our simulations, we take a Galactocentric rectangular co- 
ordinate system, where the x and y-axes are orthogonal in 
the Galactic plane and z-axis perpendicular to the plane. 
The dynamic status of a pulsar at any time can be described 
by {x,y, z;Vx,Vy,Vz), where x, y and z are position coor- 
dinates, and Vx, Vy and Vz the velocities along each axis. 
Obviously, the Galactocentric radius R = ^/x'^ +y'^. The 
distance between the Sun and the Galactic center is taken 
as Ro = 8.0 kpc. Simulation results should not depend on 
coordinate system taken. Before calculating pulsar positions 
at any time t, the initial positions and velocities of pulsars 
and the Galactic acceleration should be discussed. 



2.1 Initial height distribution 

Initial positions of pulsars should follow the positions of their 
progenitors, e.g. OB stars, or the distribution of supernova 
remnants. Maiz-Apellaniz (2001) has listed all previous de- 
terminations of scale heights of the OB star disk in either an 
exponential or a Gaussian height distribution. He obtained a 
solid measurement of the vertical distribution of B-B5 stars 
in solar neighborhood as 



f^^ini) = 



(1) 



/27r/lini 

with a scale-height of /iini = 63 pc from Hipparcos data, 
though a single-component, self-gravitating, and isothermal 
disk model with a sech function is as good as the Gaussian. 
Here wo ignore the halo component which is about 4% of 
stars, and it is not clear whether this component is related to 
MSPs. In our simulations, we assume this height distribution 
to be adequate for any places in the disk. 



We also tried an exponential distribution Pz{zini) = 
(1/2/iini) exp(— |zini|//iini) with a scale-height hini ~ 60 pc 
and a flat distribution in infinitely thin disk (aini = 0). We 
found that our results are not sensitive to different initial z- 
distributions, as all of them lead to the same or very similar 
dynamic behavior after one million year. In the following we 
will use the Gaussian distribution in Eq.(l), and assume it 
is valid at all Galactocentric radii. 



2.2 Initial radial distribution 

The radial density distribution of newly born neutron stars 
is not clear at all. Narayan (1987) made a Gaussian fit to 
the observed radial density profile obtained by Lyne, Manch- 
ester & Taylor (1985) as using a normalized function 



1 

64^ 



exp[ 



- g; J 



(kpc) 



(2) 



where pr is radial density. Bailes & Kniffen (1992) pointed 
out that another function with pH — » as i? — » is also con- 
sistent with the observed data. However, most of the sub- 
sequent work took Narayan's fit as a starting distribution 
in R (e.g. Lorimcr et al. 1993; Johnston 1994; Mukliorjee 
& Kembhavi 1997). While, this distribution reflects the ob- 
served distribution in R, rather than the initial distribution 
which most probably is related to the exponential distribu- 
tion of OB stars in the Galactic plane, as expected in our 
Galaxy and other spiral galaxies (Bahcall 1986; Paczyiiski 
1990). Note here that the radial probability distribution (Pr), 
which is used for discussions below and plots in Figure 1, 
should have the form of 

2nR pRdR 



Pr dR 



J^2nR PR dR' 



(3) 



For Narayan's function' , Pr = || exp[— (■^)'^]. In fact, this 
function has a natural deficit for _R — > which can be ev- 
idently seen from the observed data that Johnston (1994) 
tried to explain. 

In our sirrmlations, wo have tested a few probability 
functions, as shown in Figure. 1. 
(1) The exponential distribution: 



PRiRini) = 



Rc 



■exp(--^ — )' 



(4) 



hereafter i?ini ~ V^ni + iyini ^^^'^ initial Galactocentric 
radius. The re-scaled characteristic radius i?exp =4.7 kpc'^, 
as given in Hartman ct al.(1997); 
(2) The Gamma distribution: 

PR{Rini) = aR exp(-^^), (5) 

with i?cxp = 4.5 kpc and 0.4 kpc< -Rini < 25 kpc follow- 
ing Paczyiiski (1990). This radial distribution corresponds 

^ We note that in literature some authors claimed to use the 

Gaussian radial distribution of Narayan (1987), but as a mistake, 
they used the form of pr{R) to represent Pr. The prefactor of 
Pr, R/32, is R dependent, while that of Pr{R) is a constant of 
l/(647r). 

■t The authors implicitly used ilo = 8.5 kpc in context. Here we 
scaled it to ilo = 8.0 kpc. Same are the following cases 3, 4. 
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Figure 1. DifTcrcnt forms of radial distribution plotted against 
radius R. A uniform density distribution is a linearly increasing 
function Pn against R. The vertical line at -R = 0.4 kpc indicates 
the minimum R in our simulation. 



to exponential radial density. Here ur = 1.0683. See also 

Gonthier et al. (2002); 

(3) The Gaussian distribution, 



PR(i?ini) = 



2ttR, 



(6) 



with a rescaled Rg = 4.5 kpc following Lorimer et al. (1993). 
Many authors (e.g. Hartman et al. 1997) used this distribu- 
tion afterwards; 

(4) The offset Gaussian distribution: 



Pfl(ilini) = 



27rJ?, 



exp[- 



{Rini — Ros)^ 
2R^ 



(7) 



with rescaled _Rg — 1.7 kpc and RoS = 3.3 kpc given by 

Hartman et al. (1997); 

(5) The distribution from Narayan: 



PRiRini) = exp( 



(8) 



here Rg = 4.5 kpc taken from the normalized factor of the 
Gaussian distribution; 
(6) The linear distribution: 



Pfl(iii„i) = 



2Eini 



(9) 



H2 — ill 

so that in a given range of Galactocentric radii Ri < 7?ini < 
R2 (we took Ri — 0.4 kpc and R2 = 25 kpc) the pulsar 
surface density to bo a constant. 

We have simulated pulsar motions in all above cases and will 
make comparison on the evolved radial distribution, while 
for detailed studies, we will concentrate on the Gamma dis- 
tribution (Eq. 5). 



2.3 Initial velocities 

The initial velocities can be written as iTini = «birth + Crot, 
where t/birth is the birth velocity and Viot the rotation veloc- 
ity of the Galax;y. The physical origin of the birth velocity 
is not clear hitherto, which might be generated from the 
disruption of a binary (Gunn & Ostriker 1970), rocket ef- 
fects (Helfand & Tademaxu 1977) or asymmetric supernova 
explosion (Dewey & Cordes 1987). We will not go further 



into details but assume that the birth velocity is isotropic 
and follows a Maxwellian distribution 



Pv («bir 



'^birth 



exp(- 



2^birth' 



(10) 



where Ubirth is ID velocity dispersion, corresponding to the 
mean velocity of ^/%pKa■hiTt\l or a 3D velocity dispersion of 
v^Cbirth- We will show the dynamical results of cTbirth =100, 
200, 300 and 400 km in our simulations^ . 

2.4 Galactic potential and acceleration 

Pulsars accelerate in the Galactic gravitational potential, as 
the acceleration being g = — V(/>. The potential 4> is deter- 
mined by the mass distribution in our Galaxy, V'^i^ = AirGp, 
where G is the gravitational constant, and p mass density 
as a function of position. One can find that 



In fact, Eq. 11 can be rewritten as 



1 

"4^ 



dg. 1 djRgR) 
dz R OR 



(11) 



(12) 



in a cylindrical coordinate system. Here gz is the accelerar 
tion in z direction and gn in R direction. In order to sim- 
ulate pulsar dynamics, it is very important to find the best 
representation of the Galactic potential. 

2.4-1 Adopted Galactic potential model 

We used the potential model of disk/bulge originally pro- 
posed by Miyamoto & Nagai (1975), 

GMi 



ct>i{R,z) 



7?2 + 



+ («2 + 6; 



1/2 ' 



(13) 



here at, hi and Mi are model parameters, and i = 1 stands 

for the bulge and i = 2 for the disk. Using the density dis- 
tribution of the halo given by Kuijken & Gilmore (1989), 

Pc 



P = 



1 + (r-/rc)2 



(14) 



here r 



for the halo component as 



-■y^+z^, Paczyhski (1990) derived the potential 



GMe 



-I atan — 

r rc 



(15) 



where Mc = -iivpcrl and pc and Vc are model parameters. 
Paczyhski (1990) determined all the parameters (see Table 
1) for the disk, bulge and halo potentials to make a good 
agreement with the rotation curve, local volume density and 
the column density between z = ±700 pc. He took Ro = 8.0 
kpc. This potential expression has been used for simulations 
of pulsar motions by Hansen & Phinney (1997), Cordes & 
Chernoff (1998) and Gonthier et al. (2002) . Since the rota- 
tion curve for R < 0.4 kpc can not be well described by the 
formula of the potentials, as Sofue & Rubin (2001) show, all 
our simulations have been limited to R> 0.4 kpc. 

§ Note here 1 km 1 pc Myr~^. 
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Table 1. Potential parameters given by Paczynski (1990) with 
subscripts 1 and 2 corresponding to bulge and disk and c to halo. 



ai (kpc) 





bi (kpc) 


0.277 


Ml (Mq) 


1.12x10^° 


a.2 (kpc) 


3.7 


b2 (kpc) 


0.20 


M2 (Mo) 


8.07x1010 


Vc (kpc) 


6.0 


Mc (Mq) 


5.0x101° 



_ R=8kpc _ 














oCI87 .DB98 




«N090 

— 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 


xP90 aF98 




- R=4kpc 


— 1 — 1 — 1 — 1 — 1 — 1 — 1 — V 


— 1 




i-'"^^^^ 





12 3 4 5 

z (kpc) 

Figure 2. Acceleration calculated for i? = 4 kpc and 8 kpc 
for the models by Carlberg & Innanen (1987: CI87), Dehnen & 
Binney (1998: DB98), Ferriere (1998: F98), Kuijken & Gilmore 
(1989: KG89), Narayan & Ostriker (1990: NO90) and Paczynski 
(1990: P90). 



2.4.-2 Galactic potential previously used 

Previously, several models for Galactic mass distribution 
and ID acceleration have been used to simulate pulsar mo- 
tions in the ^-direction. 

1. Narayan & Ostriker (1990) modeled the mass density 
in the Galaxy as a thin layer at a = with surface density 
S and a uniform halo with density of p, and obtained the 
acceleration in z-direction as, 



z + 



2p z 



(16) 



P = 



(17) 



This model was also used by Itoh & Hiraki (1994) and tried 
by Hartman & Verbunt (1995). 

2. Kuijken & Gilmore (1989) rewrote Eq.(12) as 

4nG I dz 

for disk galaxies following Mihalas & Binney (1981). Here A 
and B are the Oort constants. After considering the double- 
exponential disks and the spherical components (halo, bulge 
and corona), they modeled the distribution of matter in the 
disk via tracer stars at high z and obtained the acceleration 
in 2-direction 



Qz = 1.04 X 10" 



1.26z 



+ 0.582 



(18) 



182 

which can be further used to determine the disk surface 
density. 

Their formula of Qz has been used by Bhattacharya 
et al. (1992) and Hartman & Verbunt (1995) to simulate 




R (kpc) 

Figure 3. The contours of Qz (upper panel) and gn (lower panel) 
calculated from potentials of Dehnen & Binney (1998: DB98), 

Paczynski (1990: P90), and Carlberg & Innancn (1987) with mod- 
ified parameters by Kuijken & Gilmore (1989: C187/KG89). The 
contour levels are at 2" pc Myr~^, with n= 0, ±1, ±2, ±3, and 
with notes near n=0. 



the pulsar motions in 2-directions for investigation of pul- 
sar magnetic field decay, and extended by Ferriere (1998) 
for ISM distribution studies. It seems to have been much 
better constrained at high z than that of Narayan & Os- 
triker (1990). See Figure 2 and discussions by Hartman & 
Verbunt (1995). 

3. The potential model for the disk/spheroid and nu- 
cleus/bulge, either directly from Carlberg & Irmanon (1987) 
or using parameters given by Kuijken & Gilmore (1989), has 
the form of 



GM 



Pb,n = 



GMb,„ 



2 -11/2' 

+ 62 + i^2 1 



+ ^2)1/2' 



(19) 



(20) 



which has been used by Lorimer et al. (1993, 1997), Hartman 
et al. (1997) and Mukherjce & Kembhavi (1997) for pulsar 
population synthesis. The model parameters a, b, hi, Pi and 
mass M for different components were constrained by the 
rotation curve data. 

4. Dehnen & Binney (1998) have constructed a set of 
mass models for our Galaxy using all available observational 
constraints. All of their models consist of the spheroidal 
bulge, halo and three disk components, namely the thin and 
thick stellar disks and the interstellar medium disk. The 
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Figure 4. The rotation curves at different z from the mass models 
of Dehnen & Binney (1998: DB98), Paczyiiski (1990: P90), and 
Carlberg & Innanen (1987) witii modified parameters by Kuijlten 
& Gilmore (1989: CI87/KG89). 



density for bulge and halo in the model is described by the 
spheroidal density distribution 



Ps = po 



(21) 



Here m = {R^ + q ^z^)^^"^, and the density of disks is given 
by 

R_ 

Rd Zd ' 

where po, ''o, <7, 7. ^ti ^d, 2d, -Rd and R^ are model 
parameters. One set of the best fitted parameters which we 
tried in this paper was derived by satisfying the total local 



Pd{R,z) 



^d g^p^ 

'2Zd Rd 



(22) 



and circular velocity Vc{Ro) 
(Model 2 in their Table 3). 



: 217 km s"^ with Ro = 8 kpc 



2.4-3 Discussions on Galactic potential models 

Obviously most models of gz (see Figure 2) are acceptable 
in the regions only around the sun and near the mid-plane, 
which can not be used to simulate pulsar motions farther 
than a few kpc from the Sun in our Galaxy. The change of 
acceleration (both and gn) with R and z (see Figure 3) 
must be considered. Previous results based on simulations 
with only g^ acceleration (e.g. Bhattacharya et al. 1992 ) 
therefore should be treated with caution. 

An accurate model for the mass distribution and accel- 
eration is needed for better understanding of pulsar motions 
in our Galaxy. The rotation curves for models of Dehnen 
& Binney (1998), Paczynski (1990), and Carlberg & Inna- 
nen (1987) with modified parameters by Kuijken & Gilmore 
(1989) are shown in Figure 4. We noticed that the potential 
given by Dehnen & Binney (1998) is almost the same with 
that by Paczyriski (1990), though newer constraints were 
used. Our simulations for both potentials produced very sim- 
ilar results. The main difference between these two models 
is the rotation curve peak at J? ~ 300 pc, which Dehnen & 
Binney (1998) did not but Paczynski (1990) did consider the 
rotation curve given by Burton & Gordon (1976) in small- J? 



part fitted from the HI data of Westerbout (1976). Although 
the peak at the small- i? is not unusual in spiral galaxies (So- 
fue et al. 1999), it is not clear whether and how the central 
bar affects the rotation in our Galaxy (Dehnen private com- 
munication). Therefore, in our simulations we only consider 
pulsar motions in the region of 0.4 kpc< R <25 kpc. 

We also noticed that the mass model of Carlberg & In- 
nanen (1987) even with modified parameters by Kuijken & 
Gilmore (1989) cannot produce the proper rotation curve 
at high z and largo R. The very slow change of the poten- 
tial near the Galactic center (Figure 3) is not reasonable. 
The accelerations given by Dehnen & Binney (1998) and 
Paczynski (1990) at small R and at high z (sec Figure 3) 
seems to be better constrained by observational data (sec 
Figure 4). Therefore the models of Dehnen & Binney (1998) 
and Paczynski (1990) are recommended here. Paczyiiski's 
potential is preferable due to its simple analytic form. 



2.5 Governing equations and calculation method 

The initial velocities and positions of the pulsars were gen- 
erated randomly according to the discussions above. The 
rotation velocity was calculated from the Galactic poten- 
tial. The motion of a pulsar is therefore only governed by 
the Newtonian kinetic equation as, 



drjt) 

dt 
dv (r) 

dt 



V (r) 
9{r), 



(23) 



where f{t), v ( f) and (r) arc 3D vectors. The equations 
above were solved numerically in three directions using the 
5th order Runge-Kutta method with adaptive control step- 
sizes [using the subroutine riqs] (Press et al. 1992). The 
initial stepsize is dt = 10"* Myr, and calculations of po- 
sition and velocity in 3D then continues until 0.1 Myr 
with adjusted stepsizes according to the required accuracy. 
The data are recorded every 0.1 Myr and the total energy 
Etot = {vl +Vy +vl)/2 + ct>{R, z) are checked. The position 
and velocity then are used for input of the next 0.1 Myr. 
The fluctuation of the total energy of each pulsar is always 
<~ 2% from to 2 Gyr, which we believe are good enough 
for simulation accuracy. 

Some examples of moving trajectories of pulsars are 
plotted in Figure 5, with various initial conditions. Obvi- 
ously, the different initial heights have very little influence 
on z tracks (Figure 5- A), while if the pulsar is located at 
different J?, trajectories can be very different (Figure 5-B). 
Certainly, different rotation velocities (Figure 5-C) or ver- 
tical velocities (Figure 5-D) will naturally cause different 
trajectories. 



3 RESULTS AND DISCUSSION 

Limited by the computation resources available, we simu- 
lated 2 X 10^ pulsars at the age t = 0, and put their initial 
positions and velocities according to the details discussed 
in Section 2. We then traced their motions up to an age of 
2 Gyr, and perform statistics on their locations every 0.1 
Myr. We first present results from Gaussian distributions of 
the initial height z, and Gamma-distributions for initial R, 
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Figure 5. The z-tracks of pulsars, corresponding initial values are marked on the right side of the plots. 



which will be called the standard initial conditions hereafter. 
We then compared those with the results from other initial 
conditions. 

3.1 Results from standard initial conditions 

3.1.1 The scale-height evolution 

The ^-heights of simulated pulsars arc binned with equal 
height interval every 0.1 Myr. We used two datarsets: the 
average height of the ith bin, z{i), and the number of pulsars 
in the ith bin, N(i). The Lovcnburg-Marquardt method was 
then employed to fit these two data-sets with given functions 
(p.678 of Press et al. 1992) to obtain the height distribution. 

For pulsars with an age less than 8 Myr, the height 
distribution can be well fitted by a single Gaussian function 
(Figure 6: t=2 Myr and 7 Myr), i.e.. 



N{i) = >lexp(- 



2hl ' 



(24) 



Here kg is the scale-height and A is amplitude. We found 
that kg increases linearly with t (Figure 7), which can be 

represented as 



kg = ha + at, 



(25) 



where ho and cr are fitting coefficients listed in Table 2. 

The relation between a and (Tbirth should be discussed. 
At small t, one could take 



Zt = Zi, 



- V^init 



(26) 



as a simplification, where zt is height at time t. To justify this 
approximation, we calculated the turn-over time of pulsars, 
which was defined as the time when pulsars change the sign 
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Figure 6. The height distribution can be well fitted by a Gaus- 
sian (t<8 Myr) or a Gaussian plus an exponential functions 
(t>8 Myr). Here Cbirth = 300 km is assumed. 



Table 2. Fitting parameters for different initial ID velocity dis- 
persion 
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Figure 8. The distribution of pulsar turn-over times. 
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Figure 9. The evolution of the scale- heights of Gaussian com- 
ponent [hg-. lower panel) and exponential component {h^: upper 
panel) . 



+ 0.4kpc<R<4kpc 

t 4kpc<R<8kpc 

o 8kpc<R<12kpc 

^ 12kpc<R<16kpc 




of Vz-, i.e. start to move towards the Galactic plane for the 
first time. It is approximately a quarter of the oscillation 
period for a pulsar moving up and down. As can be seen 
from Figure 8, most pulsars change their motion directions 
(i.e. the sign of v-^') at Myr, a time not sensitive to initial 
velocities. Therefore, for a Maxwellian initial 3D velocity 
distribution, it is natural that the distribution of height at 
small age t can have a Gaussian form, and that the scale- 
height increases with time as hg ~ (Tbirth t- We noticed that 
roughly CTbirth ~ cr, a relation that can be used to determine 
ID initial velocity dispersions. We emphasized that formula 
(26) is only an approximation of the first order. The effect 
of the deceleration of the Galactic potential did work to 
a certain extent, which made a slightly smaller, i.e. a <~ 
cTbirth (Table 2). 

As time increases, the central peak in the height distri- 
bution gets more and more prominent (Figure 6: t=10 Myr). 
A single Gaussian function cannot fit the distribution sat- 
isfactorally. We found that a Gaussian plus an exponential 
function provides a much better fit, i.e.: 



N{i) = Aexp(- 



B exp(- 



1^(^)1 



(27) 



where hg and he are scale-heights, and A and B are am- 
plitudes. The exponential component mainly models the 
smaller heights and the Gaussian component accounts for 
the larger ones. The scale-height of Gaussian component, 
hg, increases linearly with time until t ~ 40 Myr (Figure 9). 
When t >~ 40 Myr, the height distribution gets more 



z (kpc) 

Figure 10. The height distribution at t = 
The simulation used Ubirth = 300 km s~^. 



2 Gyr in 4 il-ranges. 



concentrated towards lower heights and hg decreases gradu- 
ally (Figure 9). When t >~ 200 Myr, the height distribution 
tends to stabilize, and both hg and he show no large vari- 
ations. We noticed that the larger the initial velocity dis- 
persion, the longer the stabilization time and the larger the 
resulting scale-height. 



3.1.2 Scale-heights at different radii 

We have shown earlier that the scale-height of z distribution 
of all pulsars, from 0.4 kpc< R <25 kpc. It is not clear 
whether the z distribution changes with R. 

The height distributions in four ranges of R does not 
show significant discrepancies (Figure 10). 



3.1.3 Evolution of the radial distribution 

We assume that all simulated pulsars are initially located at 
0.4 kpc< J? < 25 kpc with the maximum at R =4.5 kpc. 
During the simulation, we trace only those pulsars in the re- 
gion 0.4 kpc< R <25 kpc. As mentioned above, the Galactic 
potential near the Galactic centre may not be realistic. This 
might cause exotic pulsar motion trajectories (Carlberg & 
Innanen 1987) so that we ignore all pulsars moving close to 
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the Galaxy Figure 13. The change of the index a in a generalized fitting 
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Figure 12. The R-distribution (lower panel) and the normalized 
radial density (upper panel) at t=0 (initial) and t=2 Gyr. 



the centre {R <0.4kpc). Pulsars moving out to /? > 25 kpc 
are defined as "escaped" pulsars, since the Galactic gravi- 
tation is too wealc to bound them. The number of escaped 
pulsars increases until 200 Myr (Figure 11). The larger the 
initial velocities, the more likely the pulsars arc to escape. 
For cTbirth = 400 km s~^, more than 60% of pulsars escape 
after 100 Myr. Obviously a huge amount of pulsars, even 
after their radio emission turns off, have gone into Galactic 
halo or intergalactic space, which could be a very important 
ingredient of the dark-matter halo. 

We found that there is almost no evolution of the ra- 
dial density and radial probability distribution (Figure 12) if 
the initial radial distribution is a Gamma-distribution. The 
different velocities do not change the distribution although 
numbers of remained pulsars change a lot. 

Concerning the form of the radial distribution, there is 
a deficit towards the Galactic centre, as Johnston (1994) dis- 
cussed. A peak appears at ~ 4.5 kpc, not too different from 
the initial distribution. However, the surface density does 
not have such a deficit, in contrast to the observed pulsar 
density distribution newly determined by Lorimer (2004). 



3.2 Results for other initial conditions 

3.2.1 Alternative initial height distributions 

Two kinds of initial height distributions have been tried for 
simulations of pulsar motions: an exponential height distri- 
bution, and a flat distribution of all z\ni = 0. 

For those non-Gaussian initial height distribution, a 
Gaussian function did not always give a good fit when 
t < I Myr. So we tried another function for the height dis- 
tribution: 

iV« = ^exp{-(M^)"}, (28) 

where a = 1 corresponds to an exponential function and a = 
2 a Gaussian function. We examined a few cases of initial 
height distributions. As shown in Figure 13, for the initial 
Gaussian distribution of heights or the flat distribution with 
«ini = 0, a is always nearly 2, so that we can always flx a = 2 
in the fit. For the exponential distribution of initial heights, 
the index a increases from 1 to 2 gradually, which can be 
approximately described by a function 

a = 2-expi--— i. (29) 
^ \ 0.35 Myr J ^ ' 

When t is greater than 1 Myr, the height distributions 
can be described by a Gaussian or a Gaussian plus an ex- 
ponential function, as shown in Section 3.1.1. We therefore 
conclude that the pulsar height distribution is insensitive to 
the initial distribution after t>l Myr. 

3.2.2 Alternative initial radial distributions 

As we introduced in Section 2.1, alternative radial distribu- 
tions can be used for simulations: (1) Gamma, (2) Gaussian, 
(3) offset Gaussian, (4) exponential, (5) Narayan and (6) 
uniform density distribution. 

We first checked the evolution of height distribution, 
and found that all these radial distributions produce sim- 
ilar results, as Figure 14 shows. That is to say, the height 
distribution is insensitive to initial radial distribution. 

The evolution of radial distributions has been checked 
using both the density and radial distributions (see Fig- 
ure 15). Obviously pulsars at large R feel a smaller Galactic 
potential and hence are easy to escape after some years. 
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Figure 14. Evolution of scale-Iieights from simulations with dif- 
ferent radial distributions. 
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Figure 15. The density distribution (normalized pji) and the 
radial distribution (Pr) at t=0 yr (upper panels) and t=2 Gyr 
(lower panels) for various initial distributions. 

The final distribution seems to be more prominent at small 
R where the Galactic potential is much larger due to the 
presence of the Galactic bulge. 

Wc indeed see the fall-off at small R in both the den- 
sity distribution and radial distribution from our simulations 
with all kinds of initial radial distributions. However, com- 
pared to the observed deficit at small R recently determined 



by Lorimer (2004), the simulated deficit is much smaller in 
all cases of radial distributions. Note that what we simu- 
late simply moving neutron stars and take no account of 
their radio emission. It is likely that the beaming effect and 
the evolution of pulsar radio emission probably have to be 
considered carefully for a fair comparison. We can conclude 
from our simulation, however, that there should be a large 
number of evolved neutron stars near the Galactic center, 
which may no longer be observable as radio pulsars. 



4 APPLICATIONS 

All above simulations are made in ideal conditions for pul- 
sar motions in the Galactic potential. Previous authors have 
simulated currently observed pulsar populations, so our re- 
sults give a complementary image of pulsars moving in our 
Galaxy. As wc just mentioned, for detailed comparison be- 
tween the simulation results and observed sample, one has 
to consider beaming and the evolution of radio emission of 
pulsars. We understand that the selection effects on pulsar 
surveys are quite severe for a few reasons. For example, in 
any flux-limited survey, more luminous pulsars arc easy to 
be detected, even they are far away. Faint pulsars can be 
detected only if they are very nearby. It is not clear how 
the luminosity of pulsars evolves with age, though evidence 
available shows that young pulsars tend to be brighter. Due 
to dispersion smearing, only nearby millisecond pulsars are 
easy to be discovered, mostly at high latitudes. 

Nevertheless, at these two age extremes, i.e. young pul- 
sars and millisecond pulsars, we do not have to synthesize 
populations for comparison of the velocity distribution or 
height distribution. As wc do not expect the pulsar lumi- 
nosity to evolve much in 1 Myr, and also previous surveys 
have been mostly concentrated on the Galactic plane where 
young pulsars live, the sample of very young pulsars suffer 
much less of selections effects in the surveys. For millisecond 
pulsars, the height distribution is very stable after 200 Myr 
(see Figures 9 and 14), so that millisecond pulsars of any 
age should follow the same height distribution fitted by a 
Gaussian plus an exponential function. In the following we 
will discuss the height distribution of young pulsars and mil- 
lisecond pulsars. 

4.1 Young pulsars: Initial Velocity Dispersion 
derived from z distribution 

Most pulsar surveys have been conducted near the Galac- 
tic plane where young pulsars are predominantly found. As 
our simulations show, the scale-height is linearly increasing 
for young pulsars in the form of hg = ho + at for t < 8 
Myr, which may be used to determine the initial velocity of 
normal pulsars. 

We first take all pulsars with characteristic ages of 
1 Myr< r <8 Myr from the most updated pulsar cata- 
log (Manchester et al. 2004) ^. Pulsar distances are esti- 
mated by using the new electron density model of NE2001 
(Cordes & Lazio 2002). To make the consistency with 
Ro = 8 kpc, all distances have been scaled with a factor 

^ http://www.atnf.csiro.au/research/pulsar/catalogue/ 
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Figure 16. The scale-heights versus characteristic ages of cur- 
rently known pulsars. Solid lines represent linear fits. The num- 
bers of pulsars in the corresponding bins are indicated in the plot. 



of 8/8.5=0.94. We binned these pulsars into four groups ac- 
cording to their characteristic ages and obtained the scale- 
height hg of each group through Gaussian fitting (filled cir- 
cles in Figure 16). We then checked the relation between 
hg and ages, Eissuming that the true age t is equal to the 
characteristic age r. The resulting ID velocity dispersion is 
29 ± 12 km s~^, too small compared with previous results 
of several hundreds from proper motion measurement anal- 
yses (e. g. Lyne & Lorimer 1994; Hansen & Phinney 1997). 
This indicates strong selection effects in the assembled pul- 
sar sample, as we justify below. 

For young pulsars (t<8 Myr), z = ^ini + fzinit. If Zini <'~ 
0.06 kpc is ignored, a pulsar with z-velocity of 400 km s^^ 
will reach 0.4 kpc at 1 Myr, but 4 kpc at 10 Myr. Previ- 
ous pulsar surveys, like the Parkes multibeam survey which 
discovered about half of all known pulsars, mainly covered 
\b\ < 5°. If the average pulsar distance is about 6 kpc, all 
pulsars of |6| < 5° would have a \z\ < 0.5 kpc at or nearer 
than this average distance. Only pulsars with r < 8 Myr 
and a ^;-velocity less than v^ini w z/t = 0.5kpc/8Myr = 
63 km s~^ can be picked up in such a survey. The detected 
older pulsars must have rrmcli smaller z-vclocities, otherwise 
it have to be much farther away than a distance of 6 kpc. 

Pulsars younger than 1 Myr cannot move too far away 
even with a large initial a-velocity because of their small 
ages. As a result, the combined sample would not suffer 
much selection effects due to survey regions. This is why 
Lyne & Lorimer (1994) chose 1 Myr as the cut-off age of the 
sample. Now we try to use much more pulsars of t < 1 Myr 
to estimate 2-velocity. About 244 pulsars have been divided 
into three bins, with roughly the same number of pulsars 
each. The ID velocity dispersion in z direction derived from 
the scale-heights of z distributions of three sub-samples is 
175±56 km s~^, as shown in Figure 16. This corresponds 
to a 3D velocity dispersion of 303 ± 97 km or a mean 
velocity of 280 ± 96 km s"^ 

To compare with previous estimates of velocities, we 
should use Taylor & Cordes (1993) model for pulsar dis- 
tance, and should not include the pre-scaling factor from _Ro. 
We then found the ID velocity dispersion of 187±64 km s^^, 
or a 3D velocity dispersion of 324±110 km s^^ and the mean 
velocity of 299±102 km s~^. These values are consistent with 
estimates given by Hansen & Phinney (1997), the mean pul- 
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Figure 17. K-S test: Cumulative probability distributions of sim- 
ulated and observed samples (upper panel), and the significance 
level for different initial velocities (lower panel). 



sar velocity of 250 km s ^, or by Lyne & Lorimer (1994), 
the mean of 450 ± 90 km s~^. 



4.2 Initial Velocity Dispersion of MSPs 

MSPs have a more complex evolutionary history than nor- 
mal pulsars. Generally, MSPs are old neutron stars spun up 
by mass and momentum transferring from the companions 
(e.g. Alpar et al. (1982)). Their "birth" velocity is either 
the velocity of the binary system or its sole velocity after 
the disruption of the binary system. Another consecutive 
problem is their ages. Characteristic ages from period and 
period derivative probably do not reflect their true ages. Are 
they chronologically old? Hansen & Phinney (1997) demon- 
strated that pulsars older than lO'^ yr show the asymmetric 
drifts. Toscano et al. (1999) confirmed the asyrmiietric drifts 
of MSPs, which justify that MSPs are really dynamically old 
enough to be virialized. There is no doubt that MSPs are 
really old objects in the Milky Way. 

From our simulation, the height distribution of old pul- 
sars (e.g. t >200 Myr) is stabilized. MSPs of all ages fol- 
low the same height distribution. For the comparison of z- 
distributions of simulated sample and observed sample, we 
did not take into account selection effects for the MSP dis- 
covery. The observed sample consists of 48 MSPs from the 
latest pulsar catalog after discarding the MSPs in globular 
clusters. The simulated sample are old pulsars within 3 kpc 
from the Sun in a number of sets of simulations with ini- 
tial ID 2:-velocity dispersion from 30 km s~^ to 180 km s~^ 
(with a step of 5 km s^^). The K-S test is employed to 
check whether these two distributions are from the same 
parent distribution. We find that (Tbirth = 60 ± 10 km 
gives the largest probability. If the velocities of MSPs fol- 
low a Maxwellian distribution, the ID velocity dispersion is 
most probably 60±10 km s~^ or the mean velocity disper- 
sion most probably 96±16 km s~^, consistent with previous 
results (e.g. Lyne et al. 1998). 



Pulsar motions in our Galaxy 11 



5 CONCLUSIONS 

We presented a generalized statistical picture of how pulsars 
move in our Galaxy. The potential given by Paczyriski (1990) 
can be a good representation of mass distribution of our 
Galaxy, and has a simple analytic formula. We found that 
the final height distributions are not sensitive to the forms 
of initial z and R distributions. The height distribution can 
be well fitted by a Gaussian function within '^S Myr, and 
the scale-height increases linearly with time. After that an 
extra exponential function is required to fit the height dis- 
tribution. The height distribution gets stabilized after about 
200 Myr. 

The height distribution of pulsars younger than 
1 Myr implies directly that the mean initial velocity of 
280±96 km s"'^. Comparison of the simulated sample of 
millisecond pulsars and observed sample suggests the ID 
initial velocity dispersion of MSPs to be most probably 
60±10 km s~^, consistent with estimates given by previous 
authors. 
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